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Detailed analyses were performed for dynamic heterogeneity in supercooled liquids using three 
different numerical methods, i.e., standard molecular dynamics (MD), isoconfigurational (IC) en- 
semble analysis, and normal mode (NM) analysis. Standard MD simulations revealed that dynamic 
heterogeneity appears in supercooled liquids and becomes much more prominent as the temperature 
is lowered toward the glass transition. It was also found that the dynamical propensity obtained by 
IC analysis is coincident well with the dynamic heterogeneity determined by MD, suggesting that 
dynamic heterogeneity in supercooled liquids is strongly correlated with the initial configuration of 
the particles. To confirm this supposition, NM analysis was used to calculate the local Debye Waller 
factors (DWFs) of particles that completely disregarded the effects of momentum. Notable spatial 
correlations were observed among the dynamic heterogeneity determined by MD, the propensity 
determined by IC analysis, and DWFs determined by NM analysis. 

PACS numbers: 64.70.P-; 61.20.Lc; 61.43.Fs 



I. INTRODUCTION 

One of the long-unresolved problems in material sci- 
ence is the glass transition. In spite of the extremely 
widespread use of glass in industry, the formation pro- 
cess and dynamical properties of this material are still 
poorly understood. Numerous efforts have been devoted 
to explain the fundamental mechanisms of the slowdown 
observed in fragile glass (i.e., the sharp increase in viscos- 
ity in the vicinity of the glass transition), including clas- 
sical free volume theory and several scenarios based on 
the concept of cooperatively rearranging regions (CRR). 
However, no one has successfully identified the physical 
mechanisms behind the glass transition. Research has 
been performed for many years in an attempt to iden- 
tify events at the microscopic scale. Of the many ap- 
proaches that have been adopted, recent large-scale nu- 
merical simulations of model supercooled liquids near the 
glass transition and direct observation of colloidal parti- 
cles in colloid glasses successfully revealed the concept of 
dynamic heterogeneity, which resembles CRR. Dynamic 
heterogeneity means that the dynamic characteristics 
(i.e., particle displacements and local structural relax- 
ations) are non-uniformly distributed throughout space. 
Researchers have investigated dynamic heterogeneity in 
colloidal systems and have used simulations to exam- 
ine it in Lennard- Jones (LJ) [2, 3] and soft-core systems 
[1, [Bl, @, 0]5 dynamic heterogeneity is found in a wide 
range of glass-forming substances. It is anticipated that 
understanding of the mechanisms of dynamic heterogene- 
ity will lead to a better understanding of the slowdown 
of dynamics in those substances near the glass transition. 
However, there is currently no common and confidently 
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held description of the role of dynamic heterogeneity in 
the glass transition process. 

In 2004, Harrowell et al. proposed the use of an 
isoconfigurational (IC) ensemble in simulations to dis- 
tinguish between the static and dynamic contributions 
to the mechanisms of dynamic heterogeneity. This ap- 
proach has prompted a great deal of research based on 
the IC method, resulting in investigations of soft-core 
particle svstems 0, i, K [O, [3, LJ particle systems 
[H, [3, flHI. [lol , lattice gas systems [I^ , water molecule 
systems [18|, and other simple models of glass- forming 
materials. The IC method employs molecular dynam- 
ics (MD) in the following manner. A large number of 
samples, called the IC ensemble, are prepared with par- 
ticles in identical initial configurations. The momenta 
of the samples are then randomly selected to match the 
Maxwell Boltzmann distribution for the temperature of 
interest, and their displacements are monitored. If, for 
example, there is no correlation at all between the ini- 
tial configuration of the particles and their displacements 
from the initial configuration, it is inferred that the mean 
squares of the displacements of the particles in the en- 
semble well are approximately equal. On the other hand, 
if there is a large disparity in the mean squares of the IC 
particle displacements, resulting in a heterogeneous spa- 
tial distribution, it is possible that that heterogeneity 
was strongly influenced by the initial particle configura- 
tion. Therefore, we can consider the distribution of the 
mean square displacement of the ensemble (called the 
dynamic propensity because it indicates the freedom of 
the particles to move about) as an index of how much the 
dynamics are influenced by the initial configuration. Har- 
rowell et al. experimented with two-dimensional soft-core 
particle systems and reported that particles with differ- 
ent propensities exhibited heterogeneous spatial distri- 
butions 0. In a subsequent study, they tried to identify 
the characteristics of local static structures that deter- 
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mine the local motilities of the system. By assuming that 
the magnitude of the constraints imposed by surrounding 
particles strongly influences mobility, they investigated 
whether or not the free volume, the potential in energy- 
minimum structures, or the local Debye Waller factors 
(DWFs) could be used to predict the spatial distribution 
of mobility. It was reported that the spatial distribu- 
tion of DWFs agreed well with the spatial distribution of 
mobility [T^. 

These results indicate that as supercooling proceeds, 
static structures come to strongly influence particle mo- 
bility. This suggests that the spatial distribution of par- 
ticle mobility can be predicted from the distribution of 
DWFs, which reflects differences in the local microstruc- 
tures. Harrowell et al. employed short-term MD to cal- 
culate DWFs. Vibrational motion on the time scale de- 
fined for DWFs can be analyzed using the standard vi- 
brational analysis used in simulations of solid crystals. 
Since normal mode analysis only requires knowing the 
interactions between particles and the initial configura- 
tion of the particles, it is possible to make a complete 
prediction of DWFs while completely neglecting the in- 
fluence of the momenta of the particles. In this study, 
the propensity analysis of Harrowell et al. was extended 
to a three-dimensional system; as a result, the system 
was verified to exhibit a clear spatial heterogeneity. Nor- 
mal mode analysis was performed in an attempt to more 
clearly demonstrate the results of Harrowell et al. by 
comparing vibrational motions with propensity to mo- 
tion. This revealed a spatial correlation between active 
DWFs and propensity and, in turn, demonstrated that 
it is possible to predict dynamic heterogeneity from the 
static particle structure. 

II. MODEL AND SIMULATION METHOD 

MD calculations were carried out for a three- 
dimensional equimolar binary mixture composed of 
species 1 (small) and species 2 (large). The numbers 
of particles were Ni = N2 = 5000, and the soft-core 
potential described by Eq. ([1]) was employed in a micro- 
canonical calculation in a cube of constant volume V as 
the basic cell, surrounded by periodic boundary image 
cells. 

Vaf3{rij) = e{rij / aafB^^ , (Jafd = (cTa + Cr;3)/2, (1) 

where r^j = \ri — rj\ is the interparticle distance and 
G 1,2. The cut-off radius for the potential was set 
at 3(71. In the present paper, the following dimension- 
less units were used: length, ai; temperature, e/Zc^; and 
time, r = (micrf/e)"*^/^. The mass ratio was 1x12/ mi 
2.0. and the diameter ratio was (72/(71 = 1.2. This diam- 
eter ratio avoided crystallization of the system and en- 
sured that an amorphous supercooled state formed at low 
temperatures. The particle density was fixed at the rela- 
tively high value of p = {Ni + N2)cfI/V = 0.8. The sys- 
tem length was L = V^/^ = 23.2(7i. To reduce the com- 



putation load, a smaller system with Ni = N2 = 2500 
and L = 18.4(7i was used for IC analysis and normal 
mode analysis apart from standard MD simulations. The 
leapfrog algorithm was used with time steps of 0.005r 
when integrating the Newtonian equation of motion. 

III. RESULTS 
A. Standard molecular dynamic simulation 

The period was defined as the time when the non- 
Gaussian parameter of 0^2 (t) = 3(Ar^(t))/5(Ar^(t))^ — 
1 = 3Ari{Efii Ar,^W)/5[{Ef=i Ar|(t)>]2 - 1 shows a 
peak, providing a characteristic time at which the most 
significant dynamic heterogeneity appears. Here the 
brackets indicate ensemble averages. Arj(t) = rj{t) — 
rj{0) is the displacement of particle j during time pe- 
riod t. Figure [1] is a graphical representation of the spa- 
tial distribution of the displacement Ar{ra2) species 
1 particles during the time steps ra2' The left side of 
the figure shows the liquid state at a high temperature 
(T = 0.772), and the right side shows the supercooled 
liquid state at a low temperature (T = 0.267), which is 
around the so-called critical temperature Tc predicted by 
mode-coupling theory. Greater displacements of a parti- 
cle are indicated by more intense red color (in this study, 
the structural relaxation time was estimated at Tc ^ 0.28 
from its temperature dependence). This figure shows 
that particles with similar displacements were distributed 
in a uniform manner at high temperature. In contrast, 
in the supercooled state at the low temperature, highly 
mobile particles with high displacements and less mobile 
particles with small displacements are clustered with like 
particles, exhibiting a heterogeneous spatial distribution. 

A displacement- weighted Fourier term of the local den- 
sity Dq was defined by Eq. ([2]), and its correlation 
Sd{q) = {DqD_q) /Ni was calculated to quantify the 
dynamic heterogeneity seen in Fig. 1 [1, IH: 

where Rj = {rj{Ta^) + rj(0))/2. Figure [2] shows the 
wavenumber dependences of Sd{q) for species 1 parti- 
cles obtained at different temperatures. The correla- 
tion becomes stronger with decreasing temperature at 
low wavenumbers (long distance scales). The same trend 
was found in the correlations of species 2 particles. As 
seen in the visualization example above, this seems to 
reflect the increasing tendency for the more mobile par- 
ticles to cluster together at lower temperatures and the 
fact that this tendency is observed across long distances. 
The increasing intensity of Sd{q) observed in the low-q 
region with decreasing temperature quite closely resem- 
bles the behavior of the static structure factor of density 
fluctuations when the critical point is approached as the 
correlation length diverges and the dynamics slow. 
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B. Isoconfigurational Ensemble Analysis 

Standard MD simulations have confirmed that dy- 
namic heterogeneity increases as supercooling proceeds 
at low temperatures. However, it is difficult to conclude 
using standard MD whether meaningful relationships ex- 
ist between dynamic heterogeneity and microstructures. 
Thus, the next step was to use IC analysis, which is capa- 
ble of distinguishing the effects of momentum. One hun- 
dred distinct IC ensembles were created with randomly 
set velocities, and MD calculations were carried out. The 
method of Harrowell et al. was employed. The mean 
square displacement of the IC ensemble (Ar|(t))ic = 

X]iii°(^j(^o -\- t) — '^j{to))'^ /Niso was calculated and is 
here called the dynamic propensity of particle j. Here 
the notation (. . .)ic indicates the mean of the param- 
eter for the IC ensemble, and to is the time at which 
the velocity was reset. Figure [3] presents the spatial dis- 
tribution of the dynamical propensity at time steps of 
ra2- The left side of the figure shows the liquid state 
at a high temperature (T = 0.772), and the right side 
shows the supercooled liquid state at a low temperature 
(T = 0.289). The size and color intensity of each par- 
ticle indicate the magnitude of its displacement. Large 
numbers of particles with identical propensities at high 
temperatures showed homogeneous distributions; on the 
other hand, at low temperatures, particles were observed 
to cluster according to the level of propensity, resulting 
in a rather heterogeneous spatial distribution, analogous 
to Fig. 1. 

To quantify the dynamic heterogeneity at this propen- 
sity, these data were analyzed by the method used for 
standard MD. Propensity- weighted Fourier terms for lo- 
cal density are defined as in Eq. ([3|), and 5'pro-pro(^) = 
{DqD-q) /Ni correlations were calculated: 



E 



(Ar2(r«J) 



IC 



(Ar2(r,J)ic,: 



exp[-igf •r^-(to)], (3) 



where (. . .)ic,mean indicates the mean propensity of all of 
the particles. Figure [H presents the means of the results 
for the IC ensemble for five different initial configurations 
under both the liquid condition (T = 0.772) and the su- 
percooled liquid condition (T = 0.289). This figure shows 
that the correlation significantly increases under the su- 
percooled liquid condition at low wavenumbers. Species 
2 particles show the same tendencies in this correlation. 
This indicates that the correlation reflects dynamic het- 
erogeneity at the levels of propensity observed in the su- 
percooled state, as seen in standard MD. 

The IC analysis showed that the dynamical propensity, 
which depends only on the structure in the supercooled 
state at low temperatures, continued to show the same 
level of dynamic heterogeneity; i.e., this dependence was 
not dominated by dependency on other parameters. This 
suggests that dynamic heterogeneity is increasingly influ- 
enced by the microstructure as the temperature falls. At 



high temperatures where particle motion is especially en- 
ergetic, the dynamics of particle displacement show little 
dependence on microstructure; on the other hand, at low 
temperatures, particle motion is restrained and is much 
more subject to microstructure. It may be that dynamic 
heterogeneity is much more prominent at low tempera- 
tures due to the heterogeneity of the structure itself at 
the microscopic level. However, previous research has re- 
vealed that this microstructural heterogeneity is very dif- 
ficult to measure, not only when it is considered a static 
structure factor (i.e., a bulk parameter) but also when it 
is considered a microscopic physical quantity such as free 
volume or local potential [10, UH, [3 • 



C. Normal Mode Analysis 

The particle dynamics were analyzed by normal mode 
analysis using only the static particle configurations to 
clarify the IC analysis results, as MD is performed in the 
course of IC analysis. Unlike MD, normal mode analysis 
does not require the performance of dynamic calculations 
to determine the cooperative motion. Instead, dynamic 
matrices, which can be derived only from the static par- 
ticle configurations, provide a great deal of information 
on the vibrational motions of the particles. This analy- 
sis is rigorous for harmonic oscillations of solids. In the 
case of glassy systems, particle displacements are known 
to show a plateau regime up to roughly the a relaxation 
time. This is believed to reflect the cage effect, meaning 
that particles tend to be trapped or " caged" by surround- 
ing particles and exhibit intense vibrational motion for 
this time scale [20]. Particles are believed to show this 
vibrational behavior for a comparatively long period in 
a supercooled liquid; thus, this vibrational behavior is 
analyzed here by normal mode analysis. 

The structure in a supercooled liquid condition (T = 
0.289), which was previously observed to have promi- 
nent dynamic heterogeneity, was analyzed in this study. 
First, the initial structure used in the IC analysis was 
quenched using the steepest descent method to determine 
the local minimum potential energy structure. Next, 
the dynamical matrix V, the elements of which are 
the second derivatives of the potential energy V{r) = 
Zlili Vapirij) in mass- weighted generalized co- 

ordinates ^/mr^ was found. The eigenvalues and eigen- 
vectors of V were then determined. If the number of 
particles is then V is a 3N x 3A^ Hessian matrix, and 
its elements Vij are given by 



dridvj 



(4) 



(. . .)o indicates that the values were calculated for the 
quenched structure. Once this matrix had been simpli- 
fied to a diagonal matrix, the 3^/" eigenvalues A and the 
3A^ mutually perpendicular 3A^-dimensional eigenvectors 
a corresponding to the values of A were obtained. The 
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eigenfrequencies of the matrix were then determined from 
the relationship uj = Because the above analysis 

employs second derivatives of the potential, it is essen- 
tial to avoid truncation error in the cutoff radius of the 
potential. Therefore, a modification was applied to the 
interactions between particles in Eq. ([T]) to prevent dis- 
continuities in the interparticle potential energy and in 
the force at the cutoff radius in the process of minimiz- 
ing the energy and the normal mode analysis [2l[. Below, 
we investigate in further detail the vibrational motion of 
supercooled liquids by examining the assessment indices 
often employed in normal mode analysis. 

Density of states (DOS): The vibrational frequency uo 
found by normal mode analysis is classified into one of 
several different ranges, and the number of vibrational 
modes in each range is calculated using 



DOS(cc;) 



^ 3N 
IN ^ ^ 



3N 



(5) 



where 5 is the Kronecker delta. Figure [5] shows the fre- 
quency distribution of DOS (green lines) and of the par- 
ticipation ratio (blue dots). The participation ratio pk 
is a quantification of the extent of participation by each 
mode and is expressed by the following equation: 



N 



Pk 



(6) 



is the eigenvector of 



Here, k is the mode index, and 

particle j in mode k. "Localized mode" means that only 
a minority of the particles have very large eigenvectors. 
If all particles have eigenvectors of identical magnitude, 
Pk will be 1; if one particle has a very large eigenvector 
and if the mode is spatially localized, pk will have the 
very low value of 1/7V. Figure [5] shows participation of 
the modes and indicates low participation of the lowest- 
and highest-frequency wave modes. Typical modes in the 
frequency bands, labeled with numbers in Fig. [5l are vi- 
sualized with their eigenvectors / y/ffTJ in Fig. [6l Small 
numbers of particles at the lowest and highest frequen- 
cies had high eigenvectors, indicating the localization of 
these modes. This localization seems to be caused by the 
amorphous structure of the system. As shown below, lo- 
calized modes considerably contribute to the vibrational 
behavior of systems and are believed to have a strong 
influence on dynamic heterogeneity. 

Next, to confirm the validity of normal mode analysis, 
the dispersion relationship of this soft-core particle sys- 
tem was calculated. The dispersion relationship can be 
found from the dynamic structure factor 6'(/c,a;); how- 
ever, it can also be calculated by normal mode analysis 
using the same data. The method of Taraskin et al. was 
employed here [S^]. If an impinging plane wave with 
wavenumber k is assumed in Eq. ([7j), the contribution 
of mode j to that plane wave can be applied in Eq. (|8|) 
as a spectral density coefficient. Term A in Eq. ([7j) is 



a normalization coefficient, and a longitudinal wave and 
a transverse wave can be distinguished by whether the 
polarization vector n is parallel or perpendicular to the 
wavenumber vector /c, respectively, is the eigenvec- 
tor of particle i in mode j. This coefficient provides the 
spectral density, which is calculated using Eq. ([9|). 



Wk,i = Ancos{k ■ Vi) 



N 
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Figure [3 shows the dependences of the product of the 
spectral density and uo on wavenumber and frequency; 
i.e., it is a two-dimensional plot of the dynamic struc- 
ture factor. Deeper red color indicates a greater value. 
The wavenumber and frequency dependence of the spec- 
tral density indicated in this figure are qualitatively very 
similar to the dispersion relations in amorphous glass 
found in previous studies [l^, [l^, UBI, • The wavenum- 
ber k = Qp ^ 271 / CF = 6.28 at which the static struc- 
ture factor shows its first peak and the half-wavenumber 
k = Qp/2 ^ 3.14 are indicated by the dashed lines in the 
figure for longitudinal waves. The region between k = 
and Qp/2 is called the first Brillouin zone and is an essen- 
tial wavenumber region for investigating phonon-related 
phenomena. The figure indicates that the locations of the 
peaks in the wavenumber region are well approximated 
by a; = Csk; i.e., they exhibit a linear relationship. 

Figure 8(a) shows the dispersion relationship deter- 



mined from the wavenumber at the peak spectral den- 
sity values expressed in Eq. (|9]) . The dispersion relation- 
ship describes the characteristics of a plane wave passing 
through a material. This figure shows that the dispersion 
relationship is linear with cj = Csk in the first Brillouin 
zone, where the wavenumber is low and the wavelength 
is long. This linear relationship is characteristic of the 
dispersion relationship of phonons in a crystalline struc- 
ture; thus, it is possible to determine the speed of sound 
Cs from the gradient of this line. For comparison, a line 
representing the speed of sound found from the DSFs in 
the same type of system in a previous study is drawn 
in this figure. The relationship obtained in the present 
study is quite close to the latter line, confirming the re- 
liability of the present calculation. 

The structure examined in this study was amorphous. 
Due to its disordered structure, an impinging plane wave 
cannot be expressed by a single normal mode but rather 
by a superposition of normal mode waves with differ- 
ent frequencies. If a plane wave is introduced in an 
amorphous solid at t = 0, it is immediately decomposed 
into several different normal modes distributed within a 
narrow spectrum of frequencies for t > 0. This means 
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that the initial plane wave is dispersed by th e diso rdered 
structure and begins to attenuate. Figure [8(b)] shows 
the wavenumber dependence of the attenuation time con- 
stant Tk for the plane wave associated with wavenumber 
k, which was calculated using the half- width T{k) of the 
spectral density. The relationship r^^ ^ T{k)/2 holds. 
This figure shows that the short- wavelength vibrations of 
the high wavenumber components are damped out within 
a short time. This tendency is in reasonable qualitative 
agreement with previous findings [1^ . The time constant 
Tk is known to obey the power law, and an approximate 
relationship of r^^ ~ k'^ was obtained for this system. It 
has been suggested that the value of this exponent de- 
pends on the strength of the constraints imposed by the 
particle interactions \2§\. 

The above analysis demonstrates that the dispersion 
relationship provided by normal mode analysis agrees 
reasonably well with the results from previous research 
and supports the use of normal mode analysis to analyze 
the present system in supercooled states. Dispersion rela- 
tionships resembling those for crystals were obtained for 
long- wavelength vibrations; they mimic the behavior of 
phonons in crystals. As the more detailed analysis given 
below reveals, however, the amorphous structure of this 
system generates highly localized low- frequency modes. 
The eigenfrequencies of these modes are less than the 
minimum frequencies (longitudinal waves: 2.10; trans- 
verse waves: 0.698) of phonons corresponding to the min- 
imum wavenumber {kmin = 27r/L = 0.342) determined 
by the size L of the system, which is found from the lin- 
ear relationship between the frequency and the wavenum- 
ber. These lie outside the range of "phonons." Still, it 
is shown below that these types of low-frequency par- 
ticipating modes contribute greatly to the dynamics of 
supercooled liquids. 

Next, we attempt to quantify particle vibrational mo- 
tion to compare it with relaxation (or diffusive) mo- 
tion. The eigenmodes determine the direction and rel- 
ative magnitude of particle displacement. All particles 
oscillate at the same frequency uok in mode k. The dy- 
namics of the system can be expressed as a linear com- 
bination of all of the independent harmonic vibrational 
systems. The 3A^-dimensional eigenvector corresponding 
to eigenfrequency uJk is written as a^, and the vector cor- 
responding to particle j is written as a^. The vibrational 
amplitude vector rj^ for particle j in mode k can be writ- 
ten in the general form rjj = C^a^ / -^/rfij cos ujkt. is 
the amplitude common for all particles used to determine 
the mutual magnitude of mode k among all of the modes. 
This analysis assumes that the modes are harmonically 
independent. A thermal energy of kBT/2 is allotted 
to each mode; i.e., E^^i(^j^fe)((^i = ^bT/2. 



From this. 



we find {cos^Ukt) = \, ^jidj)'^ = 1, yield- 
mg (C^)^ = 2kBT/ujl and determining the amplitude 
for each mode. This shows that the common ampli- 
tude increases as the eigenfrequency of the mode in the 
time region in which the energy conditions are satisfied 



decreases. Considering the superposition of the eigen- 
modes for eigenfrequencies, with the exception of the 
three lowest- frequency modes that express pure trans- 
lation, the local DWF, indicating the magnitude of vi- 
bration of particle j, is calculated using: 



DWF^ 




3N-3 



J rj. 3N-3 k 



rrij 



Figure 9(a)| show s the spatial distribution of propen- 
sity, and Fig. |9(b) shows the spatial distribution of the lo- 
cal DWFs. These quantities were calculated for an iden- 
tical initial structure in a supercooled liquid state at a 
low temperature (T = 0.289). These are contour plots in 
the xy plane, as if the system box had been sliced with 
the z axis as its center; more intense red color indicates 
a greater value of the DWF. 

Both the propensity and DWF distributions indicate 
that particles with high and low values of these param- 
eters tended to cluster with like particles, as indicated 
by the non-homogeneous spacing of the particles in the 
figure. It is also noteworthy that areas of particles with 
high or low values show some notable spatial correla- 
tions between propensity and DWF. The same tendency 
was found for five different initial structures. To quan- 
tify the spatial correlation between them, Fourier compo- 
nents for local density weighted by propensity (Eq. (|3])) 
and weighted by DWF (Eq. (pT]) ) were defined, and these 
cross-correlations 5'pro-DWF(^) = {DqD'_q) /Ni were cal- 



culated: 



DWFrr 



■exp[-igr-r^-(to)], (11) 



where DWFmean is the mean value for all particles and 
Vjito) is the initial location of particle j prior to quench- 
ing. Figure [10] shows a comparison between the cross- 
correlation and the usual static structure factors (the 
green circles) for species 1. The cross correlation at the 
peak when the wavenumber Qp ^ 2it / cri = 6.28 is almost 
as large as the static structure factors, while the cross cor- 
relation is much larger than the static structure factors 
at lower wavenumber s. The same trends were found for 
species 2. This refiects a significant spatial correlation 
between the clusters of particles with high/low propen- 
sity values and the clusters of particles with high/low 
DWF values. The result that propensity and DWFs are 
spatially correlated agrees with the finding of Harrowell, 
et al. In the present study, however, it was possible 
to clearly verify the structural dependence of the propen- 
sity using the values of DWFs evaluated by normal mode 
analysis using only static particle configurations. 

Another advantage of normal mode analysis is that 
because the oscillation intensity is known for every fre- 
quency, it is possible to calculate DWFs for every fre- 
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quency band and determine the bands for which the os- 
cillation is strongly correlated with propensity. There- 
fore, the total frequency range was divided into four 
bands every 5[l/r], and the DWFs associated with those 
bands were calculated and displayed graphically. Figure 
[TTI shows the DWF for each frequency band. The range 
of colors in the contours varies with the DWF in the fre- 
quency bands; these results reveal that the DWFs in the 
low-frequency bands contribute more to the total DWF 
shown in Fig. |9(b)| A comparison of propensity and DWF 
for each frequency band shows that the calculated DWF 
for the low-frequency bands, < < 10, has a stronger 
correlation with propensity. 

IV. CONCLUSION 

In this study, we performed a standard MD simulation, 
IC ensemble analysis, and normal mode analysis to shed 
more light on the nature of dynamic heterogeneity in su- 
percooled liquids. This was done to gain a better under- 
standing of dynamic heterogeneity, which is anticipated 
to lead to an improved understanding of the fundamental 
mechanisms of glass transition phenomena. This investi- 
gation emphasized the effect of particle microstructures 
on dynamic heterogeneity. 

A calculation using standard MD revealed that dy- 
namic heterogeneity becomes much more prominent at 
low temperatures in supercooled liquids. IC analysis re- 
vealed that the dynamical propensity is coincident well 
with dynamic heterogeneity in supercooled liquids. This 
suggests that dynamic heterogeneity in supercooled liq- 
uids is strongly dependent on the initial configuration of 



the particles. Normal mode analysis that completely dis- 
regarded the effects of momentum was used to calculate 
the local DWFs of particles. A spatial correlation was 
found between the propensity determined by IC analysis 
and DWFs, and a strong dependence of dynamic het- 
erogeneity on particle configuration was confirmed. Fur- 
ther comparisons between propensity and DWF suggest 
that if static information about particle configuration can 
be obtained, DWFs can be calculated, which might al- 
low prediction of long-term relaxation dynamics to some 
extent. Analysis of the frequency bands of DWFs re- 
vealed that low-frequency modes have a stronger corre- 
lation with the spatial distribution and with propensity. 

The IC analysis and normal mode analysis were con- 
ducted using identical initial structures. MD was used 
in the IC analysis to calculate the propensity, while only 
static configurations were used to determine the DWFs 
in the normal mode analysis. These calculations revealed 
that DWF had a spatial correlation with the dynamic 
heterogeneity obtained in the IC analysis, suggesting that 
dynamic heterogeneity can be predicted from the initial 
structure. 

A recent study by Harrowell et al. found a correlation 
between irreversible reorganization of particles and the 
low-frequency soft modes in supercooled liquids consist- 
ing of a two-dimensional soft-core particle system poj . 
The irreversible reorganization employed in that study 
permits the evaluation of irreversible displacements re- 
lated to long-term structural relaxation; thus, it is de- 
fined somewhat differently from propensity. However, 
our present results coincide with those results, especially 
in the strong dependence of regions of mobility on the 
low-frequency modes. 
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FIG. 1: Spatial distributions of particle displacement Arj{Ta2) at T = 0.772 and 0.267 for time steps =1-5 and 272, 
respectively. Diameters of the depicted particles correspond to the magnitude of Arj{Ta2) plotted at the particle 

locations at t = 0. Each graph shows a full simulation system that contains 10,000 particles and has a linear dimension of 
23.1(71. 
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FIG. 2: Temperature and wavenumber dependence of Sniq) of species 1 particles. 





FIG. 3: Spatial distribution of propensity (Ar|(Ta2))ic at T = 0.772 and 0.267 for time steps = 1-75 and 143, respectively. 
Diameters of the depicted particles correspond to the magnitude of (Ar|(ra2))ic and are plotted at the particle locations at 
t — to. Each viewgraph shows a full simulation system that contains 5,000 particles and has a linear dimension of 18.4cri. 
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FIG. 4: Dynamic correlation of propensity in the liquid state and the supercooled state. Dots indicate mean values calculated 
from five different initial locations, and error bars indicate standard deviations. 
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FIG. 5: DOS and participation ratio. Green lines represent DOS frequency distributions, and blue dots represent participation 
ratios of modes. The inset on the upper right is a close-up of the participation ratio in the low-frequency region of the spectrum. 
The eigenvectors of modes labeled with numbers are plotted in Fig. [6l 
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(c) 

FIG. 6: Spatial distributions of eigenvectors at various frequencies, (a) Low- frequency localized mode, where u — 0.154 and 
Pk — 0.00946. (b) Medium- frequency delocalized mode at a; = 7.000 and pk — 0.992. (c) High-frequency localized mode, where 
ix) — 18.99 and pk — 0.00120. The origins of the vectors are located at the particle positions after quenching. There were 5,000 
particles, and the system length was 18.4cri in each case. 
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FIG. 7: Wavenumber and frequency dependences of spectral density ooSk{(^)- Intense red color indicates higher values, (a) 
Excitation by longitudinal waves. Dashed lines represent k = Qp ~ 27r/cr = 6.28 and k = Qp/2 = 3.14. (b) Excitation by 
transverse waves. 




FIG. 8: (a) Dispersion relationship under initial excitation of the first Brillouin zone by both longitudinal and transverse waves 
in a planar wave. The solid line denotes a straight line based on the speed of sound found from dynamic structure factors, (b) 
Wavenumber dependence (double log scale) of reciprocal of attenuation time constant, t^^. Solid line indicates second power 
of the wavenumber, for reference. 
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(a) (b) 

FIG. 9: Spatial distributions of propensity and DWFs calculated from identical initial structures. These are represented by 
contours of the system box in the xy plane centered on the z axis. Intense red coloring indicates higher values. The system 
contained 5,000 particles in a periodic boundary cubic box with a side length of 18.4cri. (a) Spatial distribution of propensity, 
(b) Spatial distribution of DWFs. 




FIG. 10: Comparison of cross-correlation and static structure factor (T = 0.289). Blue symbols indicate cross-correlations, 
and green symbols indicate static structural factors. Error bars indicate standard deviations of values calculated using five 
different initial locations. 
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(c)10 < < 15 (d)15 <uj <20 



FIG. 11: Spatial distributions of DWF in each frequency band calculated for identical particle configurations at T = 0.289. 
Contours of system box sliced in xy plane centered on the z axis. Intense red coloring indicates higher values. The system 
contained 5,000 particles in a periodic boundary cubic box with a side length of 18.4cri. Different ranges of colors are used in 
each DWF distribution, (a)-(d): Calculated spatial distributions for DWFs in the indicated frequency ranges: (a) < cj < 5, 
(b) 5 < cj < 10, (c) 10 < < 15, and (d) 15 < < 20. 



